/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package FEMSolver;

import BounderyPack.Boundery;
import BounderyPack.BounderyNodePair;
import BounderyPack.DirihletPair;
import FEMShapeFunctions.FEMShapeFunction1;
import HashMatrix2.HashMatrix2;
import HashMatrix2.LUFactorization;
import Integration.GaussQuadrature;
import Integration.IFunction2D;
import MathLib.triangulation.SaveStructures.Node;
import MathLib.triangulation.SaveStructures.Triangle;
import MathLib.triangulation.SaveStructures.TriangleL;
import java.io.BufferedOutputStream;
import java.io.OutputStream;
import java.io.PrintWriter;
import java.util.ArrayList;
import java.util.logging.Level;
import java.util.logging.Logger;

/**
 *
 * @author mark_o
 */
public class LagrangeMulti {
    GaussQuadrature gau = new GaussQuadrature(100);
    HashMatrix2 A;
    double[] b;
    protected ArrayList<Boundery> NapivList;
    protected ArrayList<Node> NodeListAll;
    protected ArrayList<Triangle> TrList;
    public OutputStream outPath;
    DirihletPair [] mas;
    double rho=0.005;
    double eps=0.0001;
    int Niter;

    double[] Ljambda;
    private double H=20.0;

    public LagrangeMulti(HashMatrix2 A, double[] b, ArrayList<Boundery> NapivList, ArrayList<Node> NodeListAll, ArrayList<Triangle> TrList, OutputStream outPath, DirihletPair[] mas) {
        this.A = A;
        this.b = b;
        this.NapivList = NapivList;
        this.NodeListAll = NodeListAll;
        this.TrList = TrList;
        this.outPath = outPath;
        this.mas = mas;
        Niter=0;
    }

    /*public LagrangeMulti(HashMatrix2 A, double[] b, ArrayList<Boundery> NodeList, OutputStream outPath, DirihletPair[] mas) {
    this.A = A;
    this.b = b;
    this.NodeList = NodeList;
    this.outPath = outPath;
    this.mas = mas;
    }*/


    public void calc(){

        /*for (int i = 0; i < b.length; i++) {
        for (int j = 0; j < b.length; j++) {
        System.out.print(A.getValue(i, j)+"\t");
        }
        System.out.println();
        }*/

        LUFactorization lu = new LUFactorization(A);
        //double[] u0=lu.calculate(b);
        double[] u0 = new double[b.length];

        

        Ljambda = new double[b.length];
        for (int i = 0; i < Ljambda.length; i++) {
            Ljambda[i]=1;
        }

        boolean[] gamma = new boolean[b.length];

        for (int i = 0; i < gamma.length; i++) {
            gamma[i]=false;
        }

        for (int index = 0; index < NapivList.size(); index++) {
            Boundery boundery = NapivList.get(index);

            for (int iter = 0; iter < boundery.getNodeList().size(); iter++) {
                BounderyNodePair temp = boundery.getNodeList().get(iter);
                Node vA = temp.getA(),vB = temp.getB();
                gamma[vA.Number] = true;
                gamma[vB.Number] = true;
            }
        }



        //GaussQuadrature gau = new GaussQuadrature(50);
        double[] tmpB = new double[b.length];
        while(true){
            Niter++;
            //for (int i = 0; i < TrList.size(); i++) {
            //    TriangleL tr = (TriangleL) TrList.get(i);
            //    double[] tmpNayman = getNaymanCond(tr);
            //    tmpB[tr.A-1]+= tmpNayman[0];
            //    tmpB[tr.B-1]+= tmpNayman[1];
            //    tmpB[tr.C-1]+= tmpNayman[2];
            //}
            tmpB = getNaymanCond();

            double[] newB = new double[b.length];
            for (int i = 0; i < newB.length; i++) {
                newB[i] = b[i]+tmpB[i];
            }
            double[] u1 = lu.calculate(newB);


            for (int i = 0; i < Ljambda.length; i++) {
                if (gamma[i]) {
                    Ljambda[i] = Math.max(0, Ljambda[i]-rho*(u1[i]-H));
                }else{
                    Ljambda[i] =0.0;
                }
                
            }

            if(check(u0, u1)){
                u0=u1;
                break;
            }
            u0=u1;
        }



        System.out.println(Niter+"\n==================");
        System.out.println("розвязок\n===================================================================");
        for(int i=0;i<b.length;i++){
            System.out.println((u0[i]));
        }

        //print result to file

        PrintWriter pw =  null;
        try {
            pw = new PrintWriter(new BufferedOutputStream((outPath)));
            pw.println(u0.length + mas.length);
        if(mas.length>0){
            int j;
            for (j = 0; j < mas[0].getNodeNumber(); j++) {
                pw.println(u0[j]);
            }

            for(int i=0; i<mas.length-1; i++){
                int k=0;
                pw.println(mas[i].getValue());
                while(mas[i+1].getNodeNumber() - mas[i].getNodeNumber()-k!=1){
                    pw.println(u0[j]);
                    j++;
                    k++;
                }
            }

            pw.println(mas[mas.length-1].getValue());

            for (int i = j; i < u0.length; i++) {
                pw.println(u0[i]);
            }
        }else{
            for (int i = 0; i < u0.length; i++) {
                pw.println(u0[i]);
            }
        }

        } catch (Exception ex) {
            Logger.getLogger(MyPenaltyIteration.class.getName()).log(Level.SEVERE, null, ex);
        }finally{
            pw.close();
        }

    }


    private boolean check(double[] prev, double[] next){
        for(int i=0;i<prev.length;i++){
            if(Math.abs(prev[i]-next[i])>eps){
                System.out.print("перевірка:\t");
                System.out.println(i+"\t"+Math.abs(prev[i]-next[i]));
                return false;
            }
        }
        return true;
    }

    private double[] getNaymanCond(){
        double[] res = new double[b.length];
        for (int index = 0; index < NapivList.size(); index++) {
            Boundery boundery = NapivList.get(index);

            for (int iter = 0; iter < boundery.getNodeList().size(); iter++) {
                BounderyNodePair temp = boundery.getNodeList().get(iter);
                Node vA = temp.getA(),vB = temp.getB();
                double length = ( Math.sqrt( (vA.x - vB.x)*(vA.x - vB.x) + (vA.y - vB.y)*(vA.y - vB.y) ) );
                //res[vA.Number]+=length*(H-2.0/3.0)*Ljambda[vA.Number]/2.0;
                res[vA.Number]+=length*Ljambda[vA.Number]/2.0;
                //res[vB.Number]+=length*(H-2.0/3.0)*Ljambda[vB.Number]/2.0;
                res[vB.Number]+=length*Ljambda[vB.Number]/2.0;
            }
        }

        return res;
    }

    /*private double[] getNaymanCond(TriangleL tr){
        double i = 0,
               j = 0,
               k = 0;
        double h = 0,sigma = 0;
        Node vA = NodeListAll.get(tr.A-1),
             vB = NodeListAll.get(tr.B-1),
             vC = NodeListAll.get(tr.C-1);
        boolean a1=false,a2=false,a3=false;
        for (int index = 0; index < NapivList.size(); index++) {
            Boundery boundery = NapivList.get(index);

            i=0;j=0;k=0;

            for (int iter = 0; iter < boundery.getNodeList().size(); iter++) {
                BounderyNodePair temp = boundery.getNodeList().get(iter);

                if(temp.contains(vA, vB)){
                    i=1.0;
                    j=1.0;
                    a1=true;
                }
                if(temp.contains(vA, vC)){
                    i=1.0;
                    k=1.0;
                    a2=true;
                }
                if(temp.contains(vC, vB)){
                    k=1.0;
                    j=1.0;
                    a3=true;
                }

            }
            if(i+j+k==2){
                h = Double.parseDouble(boundery.getParams().h);
                FEMShapeFunction1 fem = new FEMShapeFunction1(vA, vB, vC);
                double[] res = new double[3];
                //fem.setNodeNum(1);
                if(i!=0){
                    if(j!=0){
                        fem.setNodeNum(1);
                        res[0] = gau.GaussOnGamma(vA, vB, fem) *h;//* Ljambda[vA.Number-1];
                        fem.setNodeNum(2);
                        res[1] = gau.GaussOnGamma(vA, vB, fem) *h;//* Ljambda[vB.Number-1];
                        return res;
                    }else{
                        fem.setNodeNum(1);
                        res[0] = gau.GaussOnGamma(vA, vC, fem) *h;//* Ljambda[vA.Number-1];
                        fem.setNodeNum(3);
                        res[2] = gau.GaussOnGamma(vA, vC, fem) *h;//* Ljambda[vC.Number-1];
                        return res;
                    }
                }
                if(j!=0){
                    if(i!=0){
                        fem.setNodeNum(1);
                        res[0] = gau.GaussOnGamma(vA, vB, fem) *h;//* Ljambda[vA.Number-1];
                        fem.setNodeNum(2);
                        res[1] = gau.GaussOnGamma(vA, vB, fem) *h;//* Ljambda[vB.Number-1];
                        return res;
                    }else{
                        fem.setNodeNum(2);
                        res[1] = gau.GaussOnGamma(vB, vC, fem) *h;//* Ljambda[vB.Number-1];
                        fem.setNodeNum(3);
                        res[2] = gau.GaussOnGamma(vB, vC, fem) *h;//* Ljambda[vC.Number-1];
                        return res;
                    }
                }

            }



            if(i+j+k==2){
                h = Double.parseDouble(boundery.getParams().h);
                sigma = Double.parseDouble(boundery.getParams().g);
                double length = i*j*( Math.sqrt( (vA.x - vB.x)*(vA.x - vB.x) + (vA.y - vB.y)*(vA.y - vB.y) ) ) +
                        k*j*( Math.sqrt( (vC.x - vB.x)*(vC.x - vB.x) + (vC.y - vB.y)*(vC.y - vB.y) ) ) +
                        i*k*( Math.sqrt( (vA.x - vC.x)*(vA.x - vC.x) + (vA.y - vC.y)*(vA.y - vC.y) ) );

                //tempPyankareF = new double[]{i*h*sigma*length/2.0, j*h*sigma*length/2.0, k*h*sigma*length/2.0};
                    length*=h/6.0;
                    if(i!=0){
                        i=Ljambda[tr.A-1];
                    }
                    if(i!=0){
                        j=Ljambda[tr.B-1];
                    }
                    if(k!=0){
                        k=Ljambda[tr.C-1];
                    }

                return new double[]{
                    i/2.0*h*length , j/2.0*h*length , k/2.0*h*length
                };
            }else{
                if(i+j+k>2){
                    h = Double.parseDouble(boundery.getParams().h);
                    sigma = Double.parseDouble(boundery.getParams().g);
                    double[] res = new  double[3];
                    if(a1){
                        double length = ( Math.sqrt( (vA.x - vB.x)*(vA.x - vB.x) + (vA.y - vB.y)*(vA.y - vB.y) ) );
                        //length*;
                        res[0] = h*length/2.0;
                        res[1] = h*length/2.0;
                        res[2] = 0.0;
                    }
                    if(a2){
                        double length = ( Math.sqrt( (vA.x - vC.x)*(vA.x - vC.x) + (vA.y - vC.y)*(vA.y - vC.y) ) );
                        //length*;
                        res[0]+= h*length/2.0;
                        res[1]+= 0.0;
                        res[2]+= h*length/2.0;
                    }
                    if(a3){
                        double length = ( Math.sqrt( (vC.x - vB.x)*(vC.x - vB.x) + (vC.y - vB.y)*(vC.y - vB.y) ) );
                        //length*;
                        res[0]+= 0.0;
                        res[1]+= h*length/2.0;
                        res[2]+= h*length/2.0;
                    }
                    return res;
                }
            }
        }
        return new double[]{
                        0.0 , 0.0 , 0.0
                    };
    }
    */
}